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1  Introduction 

In  our  just  completed  research  program,  we  considered  certain  geometric  variational  meth¬ 
ods  for  problems  in  controlled  active  vision.  In  this  report,  we  will  concentrate  on  the 
work  we  performed  for  the  problems  of  registration,  warping,  and  optical  flow  as  related 
to  dynamical  tracking. 

Vision  is  a  key  sensor  modality  in  both  the  natural  and  man-made  domains.  The 
prevalence  of  biological  vision  in  even  very  simple  organisms,  indicates  its  utility  in  man-  ; 
made  machines.  More  practically,  cameras  are  in  general  rather  simple,  reliable  passive 
sensing  devices  which  are  quite  inexpensive  per  bit  of  data.  Furthermore,  vision  can  offer 
information  at  a  high  rate  with  high  resolution  with  a  wide  field  of  view  and  accuracy 
capturing  multi-spectral  information.  Finally  cameras  can  be  used  in  a  more  active  manner 
Namely,  one  can  include  motorized  lenses  mounted  on  mobile  platforms  which  can  actively 
explore  the  surroundings  and  suitably  adapt  their  sensing  capabilities. 

For  some  time  now,  the  role  of  control  theory  in  vision  has  been  recognized.  In  par¬ 
ticular,  the  branches  of  control  that  deal  with  system  uncertainty,  namely  adaptive  and 
robust,  have  been  proposed  as  essential  tools  in  coming  to  grips  with  the  problems  of  both 
biological  and  machine  vision.  These  problems  all  become  manifest  when  one  attempts 
to  use  a  visual  sensor  in  an  uncertain  environment,  and  to  fefed  back  in  some  manner  the 
information.  These  issues  constitute  a  key  thrust  in  our  research  program. 

In  this  report,  we  will  describe  in  some  detail  our  new  approach  to  the  classical  area 
of  optimal  transport.  Optimal  transport  has  appeared  in  econometrics,  fluid  dynamics, 
automatic  control,  transportation,  statistical  physics,  shape  optimization,  expert  systems, 
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and  meteorology.  In  particular,  for  the  general  visual  tracking  problem  in  controlled  active 
vision,  a  robust  and  reliable  object  and  shape  recognition  system  is  of  major  importance. 
A  key  way  to  carry  this  out  is  via  template  matching ,  which  is  the  matching  of  some  object 
to  another  within  a  given  catalogue  of  objects.  Typically,  the  match  will  not  be  exact  and 
hence  some  criterion  is  necessary  to  measure  the  “goodness  of  fit.”  The  matching  criterion 
can  also  be  considered  a  shape  metric  for  measuring  the  similarity  between  two  objects. 
In  our  work,  we  applied  concepts  from  optimal  transport  theory  to  develop  such  shape 
metrics. 

This  also  led  to  a  novel  approach  for  the  problem  of  optical  flow.  The  computation  of 
optical  flow  has  proved  to  be  an  important  tool  for  problems  arising  in  active  vision,  includ¬ 
ing  visual  tracking.  We  have  weakened  the  usual  optical  flow  constraints  with  ideas  from 
optimal  transport  and  the  theory  of  area-preserving  mappings.  This  research  represents 
our  continuing  efforts  on  the  utilization  visual  information  in  a  feedback  loop. 


2  Curvature  Flows  in  Vision  and  Image  Processing 


The  mathematical  basis  for  our  work  rests  on  two  pillars:  invariant  curvature  driven  flows 
and  geometric  variational  problems.  Here  we  will  outline  some  of  the  basic  results  on  the 
flows  which  axe  the  basis  of  the  partial  differential  equation  methods  in  controlled  active 
vision.  We  should  note  that  these  flows  themselves  are  motivated  by  ideas  from  optimal 
control  [42]. 

A  curve  may  be  regarded  as  a  trajectory  of  a  point  moving  in  the  plane.  Formally,  we 
define  a  (closed)  curve  C(-)  as  the  map  C(p)  :  S1  — *•  R2  (where  S1  denotes  the  unit  circle). 
We  assume  that  our  curves  are  have  no  self-intersections,  i.e.,  are  embedded. 

We  now  consider  plane  curves  deforming  in  time.  Let  C(p,t)  :  S1  x  [0,r)  — »  R2  denote 
a  family  of  closed  embedded  curves,  where  t  parameterizes  the  family,  and  p  parameterizes 
each  curve.  Assume  that  this  family  evolves  according  to  the  following  equation: 


f  =  «T  +  /JV 
C(p,0)  =  Co(p) 


(1) 


where  N  is  the  inward  Euclidean  unit  normal,  T  is  the  unit  tangent,  and  a  and  (3  are  the 
tangent  and  normal  components  of  the  evolution  velocity  u,  respectively.  In  fact,  it  is  easy 
to  show  that  Img[C(p,t)]  =  Img[C(m,t)],  where  C(p,  t)  and  C(w,t )  are  the  solutions  of 


Ct  =  olI  +  (3J\f  and  Ct  =  (3M, 


respectively.  (Here  Img[-]  denotes  the  image  of  the  given  parameterized  curve  in  R2.) 
Thus  the  tangential  component  affects  only  the  parametrization,  and  not  Img[-]  (which  is 
independent  of  the  parametrization  by  definition).  Therefore,  assuming  that  the  normal 
component  (3  of  v  (the  curve  evolution  velocity)  in  (1)  does  not  depend  on  the  curve 
parametrization,  we  can  consider  the  evolution  equation 

(2) 
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where  f3  =  v  ■  N. 

The  evolution  (2)  has  been  studied  extensively  for  different  functions  (5  both  from  the 
theoretical  and  applied  points  of  view.  In  particular,  it  was  introduced  into  computer 
vision  for  a  theory  of  shape  in  [39,  40,  41]. 

Two  of  the  most  important  flows  are  derived  for  (3  =  re  and  f3  —1.  In  the  former  case, 
we  will  see  that  there  is  natural  stochastic  interpretation  which  may  lead  to  an  alternative 
way  of  evolving  hypersurfaces,  which  is  one  of  the  subjects  of  our  new  research  program. 

More  precisely,  consider  the  flow 


dC  */■ 


(3) 


Equation  (3)  has  its  origins  in  physical  phenomena  [23].  It  is  called  the  geometric  heat 
equation  or  the  Euclidean  shortening  flow,  since  the  Euclidean  perimeter  shrinks  as  fast  as 
possible  when  the  curve  evolves  according  to  (3);  see  [23]. 

The  second  case  is  (3  =  1: 


dC 

dt 


=  M. 


(4) 


This  equation  simulates,  under  certain  conditions,  the  grassfire  flow  [9], 'and  is  the  basis  of 
the  morphological  scale-space  defined  by  the  disk  as  structuring  element. 

There  is  also  an  affine  analogue  of  the  Euclidean  shortening  flow  which  motivated  the 
whole  subject  of  invariant  flows.  Indeed,  in  [60],  we  show  that  the  simplest  non-trivial 
affine  invariant  flow  in  the  plane  is  given  by 


Ct  =  re^W. 


(5) 


One  can  show  that  if  C(-,0)  :  S1  — *  R2  be  a  smooth  embedded  curve  in  the  plane,  then 
there  exists  a  family  C  :  S1  x  [0,  T)  — *  M2  satisfying 

such  that  C(-,  t )  is  smooth  for  all  t  <  T ,  and  moreover  there  is  a  t0  <  T  such  that  for  all 
t  >  to,  C(-,t )  is  smooth  and  convex.  Hence  just  as  in  the  Euclidean  case,  a  non-convex 
curve  first  becomes  convex  when  evolving  according  to  (5).  After  this,  the  curve  converges 
to  an  ellipse  from  our  results  in  [60]. 

One  can  use  this  flow  to  construct  an  affine  invariant  scale-space  for  planar  shapes  [61]. 
This  in  conjunction  with  the  theory  of  differential  invariants  can  be  used  for  a  theory  of 
invariant  object  recognition. 


3  Conformal  (Geodesic)  Active  Contours 

In  this  section,  we  will  briefly  review  a  paradigm  for  snakes  or  active  contours  based  on 
principles  from  curvature  driven  flows  and  the  calculus  of  variations.  Active  contours 
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may  be  regarded  as  autonomous  processes  which  employ  image  coherence  in  order  to  track 
various  features  of  interest  over  time.  Such  deformable  contours  have  the  ability  to  conform 
to  various  object  shapes  and  motions.  Snakes  have  been  utilized  for  segmentation,  edge 
detection,  shape  modelling,  and  visual  tracking. 

In  the  classical  theory  of  snakes,  one  considers  energy  minimization  methods  where 
controlled  continuity  splines  are  allowed  to  move  under  the  influence  of  external  image 
dependent  forces,  internal  forces,  and  certain  constraints  set  by  the  user.  As  is  well-known 
there  may  be  a  number  of  problems  associated  with  this  approach  such  as  initializations, 
existence  of  multiple  minima,  and  the  selection  of  the  elasticity  parameters.  Moreover, 
natural  criteria  for  the  splitting  and  merging  of  contours  (or  for  the  treatment  of  multiple 
contours)  are  not  readily  available  in  this  framework.  In  [35],  we  propose  a  deformable 
contour  model  to  successfully  solve  such  problems  which  we  call  conformal  active  contours. 
(A  similar  approach  was  independently  formulated  in  [13,  64],  In  [13],  the  method  is  called 
geodesic  active  contours.) 

The  method  is  based  on  the  Euclidean  curve  shortening  evolution  which  defines  the 
gradient  direction  in  which  a  given  curve  is  shrinking  as  fast  as  possible  relative  to  Euclidean 
arc-length,  and  on  the  theory  of  conformal  metrics.  We  multiply  the  Euclidean  arc-length 
by  a  conformal  factor  defined  by  the  features  of  interest  which  we  want  to  extract,  and  then 
we  compute  the  corresponding  gradient  evolution  equations.  The  features  which  we  want 
to  capture  therefore  lie  at  the  bottom  of  a  potential  well  to  which  the  initial  contour  will 
flow.  Moreover,  our  model  may  be  be  extended  to  extract  3D  contours  based  on  motion 
by  mean  curvature  [35,  45]. 

Briefly,  let  g(x,  y)  be  an  image  dependent  function  tailored  to  the  type  of  feature  which 
we  want  to  capture.  For  example,  the  term  g(x,y )  may  chosen  to  be  small  near  an  edge, 
and  so  acts  to  stop  the  evolution  when  the  contour  gets  close  to  an  edge.  The  idea  now  is 
to  change  the  ordinary  arc-length  function  along  a  curve  C  =  (x(p),  y(p))T  with  parameter 
p  given  by 

ds  =  ( x\  +  yl)1/2dp, 
to 

dsg  =  (x2p  +  y2)1/2gdp. 

Then  we  want  to  compute  the  corresponding  gradient  flow  for  shortening  length  relative 
to  the  new  metric  dsg. 

Accordingly  set 

f1  dC 

L9(t) -.=  Jo  W^Wgdp. 

Then  taking  the  first  variation  of  the  modified  length  function  Lg ,  and  using  integration 
by  parts  (see  [35]),  we  get  that 

L'g(t)  =  -  -{Vg-MW)ds 
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It 
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which  means  that  the  direction  in  which  the  Lg  perimeter  is  shrinking  as  fast  as  possible 
is  given  by 

8C 

af  =  (s*-(Vj-A 0)W  (6) 

This  is  precisely  the  gradient  flow  corresponding  to  the  minimization  of  the  length  func¬ 
tional  Lj.  The  level  set  version  of  this  is 


<7> 

One  expects  that  this  evolution  should  attract  the  contour  very  quickly  to  the  feature 
which  lies  at  the  bottom  of  the  potential  well  described  by  the  gradient  flow  (7).  We  may 
also  add  a  constant  inflation  term  (see  also  [45]  for  a  more  rigorous  justification),  and  so 
derive  a  modified  model  of  (7)  given  by 


at  ~^v  IK  lv(||w|| 


)  +  v)  +  Vg-  V\F 


4  Optimal  Transport  for  Registration  and  Optical  Flow 

The  mass  transport  problem  was  first  formulated  by  Gaspar  Monge  in  1781,  and  concerned 
finding  the  optimal  way,  in  the  sense  of  minimal  transportation  cost,  of  moving  a  pile  of 
soil  from  one  site  to  another.  This  problem  was  given  a  modern  formulation  in  the  work  of 
Kantorovich  [37],  and  so  is  now  known  as  the  Monge-Kantorovich  problem.  The  problem 
of  optimal  tranpsort  has  appeared  in  econometrics,  fluid  dynamics,  automatic  control, 
transportation,  statistical  physics,  shape  optimization,  expert  systems,  and  meteorology 
[54],  It  also  naturally  fits  into  certain  problems  in  computer  vision  [25].  In  particular, 
for  the  general  visual  tracking  problem,  a  robust  and  reliable  object  and  shape  recognition 
system  is  of  major  importance.  A  key  way  to  carry  this  out  is  via  template  matching,  which 
is  the  matching  of  some  object  to  another  within  a  given  catalogue  of  objects.  Typically, 
the  match  will  not  be  exact  and  hence  some  criterion  is  necessary  to  measure  the  “goodness 
of  fit.”  For  a  description  of  various  matching  procedures,  see  [31]  and  the  references  therein. 
The  matching  criterion  can  also  be  considered  a  shape  metric  for  measuring  the  similarity 
between  two  objects. 

Registration  proceeds  in  several  steps.  First,  each  image  or  data  set  to  be  matched 
should  be  individually  calibrated,  corrected  for  imaging  distortions  and  artifacts,  and 
cleared  of  noise.  Next,  a  measure  of  similarity  between  the  data  sets  must  be  estab¬ 
lished,  so  that  one  can  quantify  how  close  an  image  is  from  another  after  transformations 
are  applied.  Such  a  measure  may  include  the  similarity  between  pixel  intensity  values,  as 
well  as  the  proximity  of  predefined  image  features  such  as  implanted  fiducials,  anatomical 
landmarks,  surface  contours,  and  ridge  lines.  Next,  the  transformation  that  maximizes  the 
similarity  between  the  transformed  images  is  found.  Often  this  transformation  is  given  as 
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the  solution  of  an  optimization  problem  where  the  transformations  to  be  considered  are 
constrained  to  be  of  a  predetermined  class.  Finally,  once  an  optimal  transformation  is 
obtained,  it  is  used  to  fuse  the  image  data  sets. 

As  we  will  explicate  below,  the  method  we  propose  for  registration  in  the  context  of 
tracking  is  based  on  an  optimization  problem  built  around  the  L2  Kantorovich-Wasserstein 
distance  taken  as  the  similarity  measure.  The  constraint  that  we  have  put  on  the  transfor¬ 
mations  considered  is  that  they  obey  a  mass  preservation  property.  Thus,  we  match  mass 
densities  in  this  method,  which  may  be  thought  of  as  weighted  areas  in  2D  or  weighted 
volumes  in  3D. 

We  will  also  describe  the  use  of  such  ideas  for  the  problem  of  optical  flow.  The  com¬ 
putation  of  optical  flow  has  proved  to  be  an  important  tool  for  problems  arising  in  active 
vision,  including  visual  tracking.  The  optical  flow  field  is  defined  as  the  velocity  vector 
field  of  apparent  motion  of  brightness  patterns  in  a  sequence  of  images.  We  have  explored 
various  constrained  optimization  approaches  for  the  purpose  of  accurately  computing  op¬ 
tical  flow.  In  our  work,  we  formulated  approaches  based  upon  the  theory  of  optimal  mass 
transport  and  area-preserving  mappings. 

4.1  Formulation  of  the  Problem 

We  now  give  a  modern  formulation  of  the  Monge-Kantorovich  problem.  Let  flo  and  fii  be 
two  subdomains  of  Rd,  with  smooth  boundaries,  each  with  a  positive  density  function,  /x0 
and  /q,  respectively.  We  assume 


so  that  the  same  total  mass  is  associated  with  fio  and  Oi .  We  consider  diffeomorphisms  u 
from  (fi0,Mo)  to  (fli,y/i)  which  map  one  density  to  the  other  in  the  sense  that 


Ho  =  \Du\  Hi  °  u,  (9) 

which  we  will  call  the  mass  preservation  (MP)  property,  and  write  u  G  MP.  Equation  (9)  is 
called  the  Jacobian  equation.  Here  \Du\  denotes  the  determinant  of  the  Jacobian  map  Du. 
In  particular,  Equation  (9)  implies,  for  example,  that  if  a  small  region  in  Q0  is  mapped  to 
a  larger  region  in  fti,  then  there  must  be  a  corresponding  decrease  in  density  in  order  for 
the  mass  to  be  preserved.  A  mapping  u  that  satisfies  this  property  may  thus  be  thought 
of  as  defining  a  redistribution  of  a  mass  of  material  from  one  distribution  Ho  to  another 
distribution  hi- 

There  may  be  many  such  mappings,  and  we  want  to  pick  out  an  optimal  one  in  some 
sense.  Accordingly,  we  define  the  Lp  Kantorovich-Wasserstein  metric  as  follows: 

dp(Ho,Hi)P  ■=  _  in,f  /  \\u{x)  -  x\\pp,0(x)  dx.  (10) 

u  €  MP  J 
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An  optimal  MP  map,  when  it  exists,  is  one  which  minimizes  this  integral.  This  functional 
is  seen  to  place  a  penalty  on  the  distance  the  map  u  moves  each  bit  of  material,  weighted 
by  the  material’s  mass. 

The  case  p  =  2  has  been  extensively  studied.  The  L2  Monge-Kantorovich  problem 
has  been  studied  in  statistics,  functional  analysis,  and  the  atmospheric  sciences;  see  [17,  8] 
and  the  references  therein.  A  fundamental  theoretical  result  [43,  10,  27],  is  that  there  is 
a  unique  optimal  u  G  MP  transporting  po  to  /h,  and  that  this  u  is  characterized  as  the 
gradient  of  a  convex  function  w,  i.e.,  u  =  Vw.  Note  that  from  Equation  (9),  we  have  that 
w  satisfies  the  Monge-Ampere  equation 

\Hw\  pi  o  (Vw)  =  pa, 

where  \Hw\  denotes  the  determinant  of  the  Hessian  Hw  of  w. 

Hence,  the  Kantorovich-Wasserstein  metric  defines  the  distance  between  two  mass  den¬ 
sities,  by  computing  the  cheapest  way  to  transport  the  mass  from  one  domain  to  the  other 
with  respect  to  the  functional  given  in  (10),  the  optimal  transport  map  in  the  p  =  2  case 
being  the  gradient  of  a  certain  function.  The  novelty  of  this  result  is  that  like  the  Riemann 
mapping  theorem  in  the  plane,  the  procedure  singles  out  a  particular  map  with  preferred 
geometry. 

4.2  Monge-Kantorovich  and  Optimal  Transport 

In  our  research,  we  focused  on  the  uses  of  ideas  from  optimal  transport  for  problems  in 
controlled  active  vision  and  visual  tracking.  However,  given  the  potential  power  of  these 
ideas  in  systems  and  control,  we  would  like  to  list  some  key  uses  of  Monge-Kantorovich: 

1.  Lyapunov  theory  is  essential  is  studying  nonlinear  system  stability  and  controller 
synthesis.  In  some  very  interesting  work,  Rantzer  [55]  has  formulated  a  dual  to 
Lyapunov’s  second  theorem.  The  idea  is  that  the  Lyapunov  function  is  regarded 
as  the  “cost  to  go”  in  an  optimal  transport  problem  and  is  dual  (in  the  sense  of 
linear  programming)  to  the  density  function  typically  studied  in  Monge-Kantorovich 
theory.  These  ideas  give  a  powerful  new  tool  in  studying  nonlinear  system  analysis. 

2.  Shape  optimization  is  another  area  of  use  for  optimal  transport.  For  example,  given 
two  densities  and  an  insulating  medium  into  which  we  place  a  fixed  amount  of  con¬ 
ducting  material  one  can  consider  the  problem  of  the  optimal  placement  of  the  con¬ 
ducting  material  to  minimize  the  heating  induced  by  the  flow.  This  can  be  put 
into  the  Monge-Kantorovich  optimal  transport  framework.  Similar  remarks  apply 
to  problems  in  compression  molding,  where  one  considers  an  incompressible  plastic 
material  being  pressed  being  two  plates  in  which  one  wants  to  track  the  air-plastic 
interface. 

3.  One  of  the  most  beautiful  uses  of  optimal  transport  is  in  meteorology,  in  particular 
semigeostrophic  models.  These  are  concerned  with  with  large  scale  stratified  flows 
with  front  formation  [17].  The  idea  is  that  meteorologists  want  to  model  how  fronts 
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arise  in  large-scale  weather  patterns.  Tracking  such  fronts  is  a  key  goal,  and  semi- 
geostrophic  equations  seem  to  give  a  reasonable  mathematical  model  for  the  creation 
of  such  fronts.  This  leads  naturally  to  optimal  mass  transport  equations. 

4.3  Background  on  Algorithms  for  Computing  The  Transport 
Map 

There  have  been  a  number  of  algorithms  considered  for  computing  an  optimal  transport 
map.  For  example,  methods  have  been  proposed  based  on  linear  programming  [54],  and 
oh  Lagrangian  mechanics  closely  related  to  ideas  from  the  study  of  fluid  dynamics  [8] .  An 
interesting  geometric  method  has  been  formulated  by  Cullen  and  Purser  [17] . 

One  very  common  method  is  to  reduce  the  L2  optimal  transport  to  a  linear  programming 
problem.  Thus  one  can  approximate  the  the  densities  po  and  po  as  weighted  sums  delta 
functions,  such  as 

Ho(x)  =  1/N  ^2  6(x  ~  xi)>  A*iO?)  =  1/JV  J2  S(x  ~  W)> 

i=l,N  i=l,N 


for  2 A  give  points  Xi, . . .  ,  xjy,  yi, ...  ,  £  Rd.  The  L2  Kantorovich  distance  is  then 

N 

Vi)2  =  infy'|yi -x^)!2,  (11) 

cr  ■ 

i 

where  the  infimum  is  over  all  permutations  on  N  letters  <x. 

This  problem  can  be  solved  as  linear  programming  problem,  by  noting  that  Equa¬ 
tion  (11)  can  be  expressed  as 

N 

inf  ^  '  CijPij 

p  : .. 

»j=i 

where  Cij  —  \xi  —  yj |2,  and  p  denotes  any  N  x  N  matrix  with  non-negative  entries,  such 
that  the  sum  of  all  columns  and  rows  equals  1  (these  are  the  “doubly-stochastic  matrices”). 
There  are  optimal  algorithms  for  general  cost  matrices  c^-,  but  to  the  best  of  our  knowledge 
there  are  no  known  optimal  algorithms  for  the  special  case  in  which  c^-  =  | X{ —yj\2-  Finally, 
note  that  even  in  the  2D  case,  typical  image  sizes  run  to  512  x  512,  and  so  the  linear 
programming  problem  can  get  to  be  quite  unwieldy. 

A  more  effective  algorithm  based  on  ideas  from  continuum  mechanics  was  proposed  in 
[8].  This  based  on  ideas  from  Lagrangian  mechanics  and  a  certain  relaxation  method.  This 
has  influenced  our  approach  discussed  below.  In  our  case  however  for  image  tracking,  we 
will  argue  that  the  most  effective  method  should  be  based  on  gradient  descent  and  the 
concept  of  “polar  factorization.”  We  discuss  this  in  the  next  section. 
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4.4  Variational  Algorithms  for  Optimal  Transport 

In  this  section,  we  describe  a  natural  solution  to  L2  Monge-Kantorovich  based  on  the 
equivalent  problem  of  polar  factorization ;  see  [10,  26,  48]  and  the  references  therein.  The 
mathematical  basis  for  the  approach  described  here  may  be  found  in  [5],  and  applications 
to  visual  tracking  in  [30].  We  will  work  with  the  general  case  of  subdomains  in  Rd,  and 
point  out  some  simplifications  that  are  possible  for  the  R2  case. 

As  above,  let  fio,  Hi  C  Rd  be  subdomains  with  smooth  boundaries,  with  corresponding 
positive  density  functions  po  and  p\  satisfying  fQo  po  —  Jq1  Mi-  Let  u  :  (fl0,  Mo)  i^-uMi) 
be  an  initial  mapping  with  the  mass  preserving  (MP)  property.  Then  according  to  the 
generalized  results  of  [10,  26],  one  can  write 

), 

u  —  (V«j)  os,  (12) 

t 

* 

where  w  is  a  convex  function  and  s  is  an  MP  mapping  s  :  (£10,  pf)  — »  (fl0,p0)-  This  is 
the  polar  factorization  of  u  with  respect  to  po-  In  [26],  just  the  case  of  area  preservation  is 
considered,  i.e.,  po  is  assumed  constant,  but  the  general  case  goes  through  as  well. 

Our  goal  is  to  find  the  polar  factorization  of  the  MP  mapping  it,  according  to  the 
following  strategy.  We  consider  the  family  of  MP  mappings  of  the  form  u  =  uo  s~l  ass 
varies  over  MP  mappings  from  (S70,  /xo)  to  itself.  If  we  consider  u  as  a  vector  field,  we  can 
always  find  a  function  w  and  another  vector  field  %,  with  div(x)  =  0,  such  that 


u  =  Ww  +  x, 


i.e.,  we  can  decompose  u  into  the  sum  of  a  curl-free  and  divergence-free  vector  field  [66]. 
Thus,  what  we  try  to  do  is  find  a  mapping  s  which  will  yield  a  u  without  any  curl,  that  is, 
such  that  u  =  Vw.  Once  such  an  s  is  found,  we  will  have  u  =  «os  =  (Vio)  o s  and  so  we 
will  have  found  the  polar  factorization  (12)  of  our  given  function  u. 

Now,  here  is  the  key  point.  As  we  discussed  above,  the  unique  optimal  solution  of 
the  L 2  Monge-Kantorovich  problem  has  the  form  u  =  Vw,  and  so  the  problem  of  finding 
the  polar  factorization  of  u  and  finding  the  optimal  Monge-Kantorovich  mapping  u  are 
equivalent.  In  essence,  our  proposed  approach  to  solve  the  Monge-Kantorovich  problem  is 
to  create  a  “rearrangement”  of  an  initial  vector  field  u  using  a  map  s,  so  that  the  resulting 
vector  field  u  =  u  o  s-1  has  no  curl. 

Finding  an  Initial  Mapping: 

We  will  describe  an  explicit  algorithm  to  solye  the  Monge-Kantorovich  problem.  So  we 
want  to  minimize  the  L2  Kantorovich-Wasserstein  distance  functional  over  MP  functions 
from  (n0)  Po)  to  (fli,  pi).  We  will  try  to  do  this  by  finding  an  initial  MP  mapping  u  and  then 
minimizing  over  u  =  nos-1  by  varying  s  over  MP  mappings  from  to  flo,  starting  with  s 
equal  to  the  identity  map.  Our  first  task  is  to  find  and  initial  MP  mapping  u.  This  can  be 
done  for  general  domains  using  a  method  of  Moser  [51,  18],  or  for  simpler  domains  using 
the  following  algorithm.  For  simplicity,  we  work  in  R2  and  assume  fl0  =  fir  =  [0,  l]2,  the 
generalization  to  higher  dimensions  being  straightforward.  We  define  a  function  a  =  a(x) 
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rx  pi 


by  the  equation 

raix)  rl  rx  r  l 

/  /  Ml  (V,y)  dy  dr)  =  I  /  fJo(y,y)  dy  drj 

Jo  Jo  Jo  Jo 

which  gives  by  differentiation  with  respect  to  x 

o!{x)  f  fj,i(a(x),y)  dy  =  f  n0(x,y)  dy. 

Jo  Jo 

We  may  now  define  a  function  b  —  6(x,  y )  by  the  equation 

rb(x,y)  ry 

a'(x)  ^!(a(x),/9)  dp  =  /  Mo (x,p)  dp, 

Jo  Jo 

and  set  u(x,y )  =  (a(x),b(x,y)).  Since  ay  =  0,  \Du\  =  axby,  and  differentiating  (15)  with 
respect  to  y  we  find 


(13) 


(14) 


(15) 


a'(x)bv(x,y)  m(a(x),b(x,y))  =  Mo (*,!/) 

\Du\niou  —  fi o, 

which  is  the  MP  property  we  need.  This  process  can  be  interpreted  as  the  solution  of  a 
one-dimensional  Monge-Kantorovich  problem  in  the  x  direction  followed  by  the  solution 
of  a  family  of  one-dimensional  Monge-Kantorovich  problems  in  the  y  direction. 

The  map  created  above  can  be  quite  crude,  and  so  one  can  use  other  methods  for 
writing  down  an  initial  mapping  including  using  an  approach  based  on  volume-preserving 
diffeomorphisms  as  suggested  in  [51]. 

Removing  the  Curl: 

Once  an  initial  MP  u  is  found,  we  need  to  apply  the  process  which  will  remove  its  curl.  We 
first  note  that  the  composition  of  two  mass  preserving  (MP)  mappings  is  an  MP  mapping, 
and  the  inverse  of  an  MP  mapping  is  an  MP  mapping.  Thus,  since  u  is  an  MP  mapping, 
we  have  that  u  —  uo  s-1  is  an  MP  mapping  if  and  only  if  s  is,  that  is,  if  and  only  if 

p0  =  |-Ds|  Mo  °  s. 

In  particular,  when  po  is  constant,  this  equation  requires  that  s  be  area  or  volume  preserv¬ 
ing. 

Next,  rather  than  working  with  s  directly,  we  solve  the  polar  factorization  problem 
via  gradient  descent.  Accordingly,  we  will  assume  that  s  is  a  function  of  time,  and  then 
determine  what  st  should  be  to  decrease  the  L2  Monge-Kantorovich  functional.  This  will 
give  us  an  evolution  equation  for  s  and  in  turn  an  equation  for  ut  as  well,  the  latter  being 
the  most  important  for  implementation.  By  differentiating  u  o  s  =  u  with  respect  to  time, 
we  get 


ut  =  ~Du  st 
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where  we’ve  abused  notation  to  define  st  :=  st°s  l.  We  need  to  make  sure  that  s  maintains 
its  MP  property.  Differentiating  Ho  =  \Ds\  //0  °  s  with  respect  to  time,  we  derive 

div(p0  st)  =  0, 

from  which  we  see  that  st,  st  and  ut  should  have  the  following  forms: 

*  =  -7-C,  (16) 

Ha 


St  =  UcJos’ 

1  , 

Ut  = - DuQ, 

Ho 


for  some  vector  field  (  on  Do,  with  div(C)  =  0  and  ((,  n)  =  0  on  dfi0,  n  being  the  normal 
to  the  boundary  of  Do-  This  last  condition  ensures  that  s  remains  a  mapping  from  Do  to 

itself,  by  preventing  the  flow  of  s,  given  by  st  =  °  s,  from  crossing  the  boundary  of 

Do-  This  also  means  that  the  range  ofu  —  uo  s-1  is  always  u( Do)  =  Di. 

Consider  now  the  problem  of  minimizing  the  Monge-Kantorovich  functional: 


J \\u-x\\2no 

J \\u\\2Ho~2  J  (u,x)ho  +  J  ||a;||2/xo- 


The  last  term  is  obviously  independent  of  time.  Interestingly,  so  is  the  first, 


J \\u\\2Ho  =  J\\uoS  '\\2Ho 

=  J\\uos-r\Ds-^os-' 

=  J IMlVo 


\u\\2Ho 


where  Ho  =  |  Ds  1 1  po  °  s  1  since  s  1  is  an  MP  map. 

Turning  now  to  the  middle  term,  we  do  a  similar  trick, 


J(u,x)h o  =  J  (uos  l,sos  x)  no 

=  J  (u  o  s-1,  s  o  s_1)  |Z?s-1|  no  °  s_1 
=  J  (u,  s)  ho, 


11 


and  taking  st  = 


o  s,  we  compute 


/t o 

\Ds\/j,0  os 


Now  decomposing  u  as  u  =  Vw  +  x,  we  have 


=  J  (Vw  +  X)  0 
=  J  (Vw,0  +  J  (x,  0 
=  J  (div(ioC)  -  wdiv(C))  +  J  ( X ,  C) 

=  /  w  (C, »)  +  f  (X,  0 

=  J  (x,  C>» 

where  we’ve  used  the  divergence  theorem,  div(()  =  0,  and  (£,n)  =  0  on  cAfl0-  Thus,  in 
order  to  decrease  M,  we  can  take  C  —  X  with  corresponding  formulas  (16)-(18)  for  st,  s4, 
and  ut,  provided  that  we  have  div(x)  =  0  and  (x,  n)  =  0  on  <9fl0.  Thus  it  remains  to  show 
that  we  can  decompose  u  as  u  =  Vw  +  x  for  such  a  x- 

Gradient  Descent:  Rd: 

We  let  w  be  a  solution  of  the  Neumann-type  boundary  problem 


,  div(ii)  =  Aw  (21) 

(Vw,  n)  =  ( u ,  n)  on  <9flo,  (22) 

and  set  x  =  ft  —  Vw.  It  is  then  easily  seen  that  x  satisfies  the  necessary  requirements. 
Thus,  by  (18),  we  have  the  following  evolution  equation  for  u: 

ut  =  — —  Du  ( u  —  VA_1div(n)) .  (23) 

At  o 
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This  is  a  first  order  non-local  scheme  for  ut  if  we  count  A-1  as  minus  2  derivatives.  Note 
that  this  flow  is  consistent  with  respect  to  the  Monge-Kantorovich  theory  in  the  following 
sense.  If  u  is  optimal,  then  it  is  given  as  u  =  Vw,  in  which  case  u  —  VA-1div(u)  = 
Vw  —  VA-1div(Vu>)  =  0  so  that  by  (23),  ut  =  0. 

Gradient  Descent:  IR2: 

The  situation  is  somewhat  simpler  in  the  R2  case,  due  to  the  fact  that  a  divergence  free 
vector  field  x  can  in  general  be  written  as  x  =  Vx/i  for  some  scalar  function  h,  where  _L 
represents  rotation  by  90  deg,  so  that  Vx/i  =  (— hy,  hx).  In  this  case,  (21)  becomes 

~\Mt  =  I  <VX/,  Vxfr)  =  J  (V/,Vh)  (24) 

where  the  decomposition  of  u  is  u  =  Vw  +  Vx/,  and  we  can  take  h  —  f.  The  function  / 
can  be  found  by  solving  the  Dirichlet-type  boundary  problem 

— div(ftx)  =  A/,  (25) 

/  =  Oon5fi0,  (26) 

which  gives  us  the  evolution  equation 

ut  =  —  DuVxA_1div(ux).  (27) 

Mo 

We  may  also  derive  a  second  order  local  evolution  equation  for  u  by  using  the  divergence 
theorem  with  (24)  to  get 


Ut  =  — —  DuVx  div(itx),  (28) 

Mo 

where  appropriate  handling  of  the  evolution  at  the  boundary  is  necessary.  If  the  boundary 
is  square,  then  one  natural  thing  to  do  would  be  to  assume  that  the  displacement  map 
u  —  x  is  periodic. 

Defining  the  Warping  Map: 

Typically  in  elastic  registration,  one  wants  to  see  an  explicit  warping  which  smoothly 
deforms  one  image  into  the  other.  This  can  easily  be  done  using  the  solution  of  the 
Monge-Kantorovich  problem.  Thus,  we  assume  now  that  we  have  applied  our  gradient 
descent  process  as  described  above  and  that  it  has  converged  to  the  Monge-Kantorovich 
mapping  uMK. 

Following  the  work  of  Benamou  and  Brenier,  [8] ,  (see  also  [27]) ,  we  consider  the  following 
related  problem: 

inf 


It 


m(£,  x)  ll'y(t)  x)  ||2  dt  dx 
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over  all  time  varying  densities  //  and  velocity  fields  v  satisfying 

^  +  div(nv)  =  0,  (29) 

A*(0,  •)  =  A*o,  ~  M1*  0  =  (30) 

It  is  shown  in  [8]  that  this  infimum  is  attained  for  some  nmin  and  vmin,  and  that  it  is 
equal  to  the  L2  Kantorovich-Wasserstein  distance  between  fj, o  and  n\.  Further,  the  flow 
X  —  X(x,t)  corresponding  to  the  minimizing  velocity  field  vjnin  via 

X(Xy  0)  =  X)  Xf  =  vmin  o  X 

is  given  simply  as 

X(x,t)  =  x  +  t  (uMk(x) -x).  (31) 

Note  that  when  t  =  0,  X  is  the  identity  map  and  when  t  =  1,  it  is  the  solution  umk  to 
the  Monge-Kantorovich  problem.  This  analysis  provides  appropriate  justification  for  using 
(31)  to  define  our  continuous  warping  map  X  between  the  densities  Hq  and  /iX. 

4.5  More  General  Cost  Functions 

We  can  also  consider  more  general  cost  functions  which  are  useful  in  image  interpolation 
and  optical  flow;  see  Section  5  below.  Indeed,  the  theory  we  described  above  extends 
formally  quite  easily  to  general  cost  functions  [5] .  Indeed,  suppose  we  are  trying  to  minimize 
a  functional 

M  —  J  $(u(x)  —  x)fio(x)  dx  (32) 

where  $  :  Rd  — >  M  is  a  positive  C 1  cost  function. 

From  the  above  argument,  we  know  that  the  evolution  of  u  :  — >  M'*  should  be  given 

as 

ut  =  -—Du-  C,  (33) 

fJ'O 

where  £  is  a  divergence  free  vector  field.  Taking  the  first  variation,  one  can  derive  an 
expression  of  the  form 

Mt  =  - J  (&(u(x)-x)X).  (34) 

Next  decomposing 


&(u(x)  —  x)  =  Vw  +  x,  (35) 
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we  have 


Mt  =  ~  J  (x  +  Vw,C) 

=  -  J  (x,C), 

by  the  divergence  theorem.  Thus  in  order  to  decrease  M  we  can  take  C  =  X>  provided  we 
show  that  div(x)  =  0  and  (x,  n)  =  0.  So  we  need  to  solve 

div($,(u(a;)  —  a;))  =  Vw, 

(Vw,n)  =  ($'(u  -  x),n)  on  dfl0 

for  w.  This  leads  to  the  following  nonlocal  equation: 

ut  =  -—Du  ■  (J  -  VA “1V-)^/(ti  -  x). 

A  similar  argument  to  that  given  above  also  will  give  a  local  gradient  descent  scheme 
in  this  case  as  well  [5]. 

5  Area- Preserving  Diffeomorphisms  and  Optical  Flow 

The  computation  of  optical  flow  has  proved  to  be  an  important  tool  for  problems  arising 
in  active  vision.  The  optical  flow  field  is  the  velocity  vector  field  of  apparent  motion  of 
brightness  patterns  in  a  sequence  of  images  [34].  One  assumes  that  the  motion  of  the 
brightness  patterns  is  the  result  of  relative  motion,  large  enough  to  register  a  change  in 
the  spatial  distribution  of  intensities  on  the  images.  Thus,  relative  motion  between  an 
object  and  a  camera  can  give  rise  to  optical  flow.  Similarly,  relative  motion  among  objects 
in  a  scene  being  imaged  by  a  static  camera  can  give  rise  to  optical  flow.  The  problem  of 
computing  optical  flow  is  ill-posed  and  so  well-posedness  has  to  be  imposed  by  assuming 
suitable  a  priori  knowledge.  For  example,  a  number  of  researchers  have  considered  a 
variational  formulation  for  imposing  such  a  priori  knowledge. 

One  constraint  which  has  often  been  used  in  the  literature  is  the  “optical  flow  con¬ 
straint”  (OFC).  The  OFC  is  a  result  of  the  simplifying  assumption  of  constancy  of  the 
intensity,  I  =  I(x,  y,  t),  at  any  point  in  the  image  [34].  It  can  be  expressed  as  the  following 
linear  equation  in  the  unknown  variables  u  and  v 


IxU  +  IyV  +  It  ■ — 0,  (36) 

where  Ix,  Iy  and  It  are  the  intensity  gradients  in  the  x ,  y,  and  the  temporal  directions 
respectively,  and  u  and  v  are  the  x  and  y  velocity  components  of  the  apparent  motion  of 
brightness  patterns  in  the  images,  respectively.  It  has  been  shown  that  the  OFC  holds 
provided  the  scene  has  Lambertian  surfaces  and  is  illuminated  by  either  a  uniform  or  an 
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isotropic  light  source,  the  3-D  motion  is  translational,  the  optical  system  is  calibrated  and 
the  patterns  in  the  scene  are  locally  rigid. 

It  is  not  difficult  to  see  from  equation  (36)  that  computation  of  optical  flow  is  unique 
only  up  to  computation  of  the  flow  along  the  intensity  gradient  V/  =  (Ix,  Iy)T  at  a  point 
in  the  image  [34].  (The  superscript  T  denotes  “transpose.”)  This  is  the  celebrated  aperture 
problem.  One  way  of  treating  the  aperture  problem  is  through  the  use  of  regularization  in 
computation  of  optical  flow,  and  consequently  the  choice  of  an  appropriate  constraint.  A 
natural  choice  for  such  a  constraint  is  the  imposition  of  some  measure  of  consistency  on 
the  flow  vectors  situated  close  to  one  another  on  the  image. 

In  their  pioneering  work,  Horn  and  Schunk  [34]  use  a  quadratic  smoothness  constraint. 
The  immediate  difficulty  with  this  method  is  that  at  the  object  boundaries,  where  it 
is  natural  to  expect  discontinuities  in  the  flow,  such  a  smoothness  constraint  will  have 
difficulty  capturing  the  optical  flow.  For  instance,  in  the  case  of  a  quadratic  constraint  in 
the  form  of  the  square  of  the  norm  of  the  gradient  of  the  optical  flow  field  [34],  the  Euler- 
Lagrange  (partial)  differential  equations  for  the  velocity  components  turn  out  to  be  linear 
elliptic.  The  corresponding  parabolic  equations  therefore  have  a  linear  diffusive  nature, 
and  tend  to  blur  the  edges  of  a  given  image.  In  the  past,  work  has  been  done  to  try 
to  suppress  such  a  constraint  in  directions  orthogonal  to  the  occluding  boundaries  in  an 
effort  to  capture  discontinuities  in  image  intensities  that  arise  on  the  edges.  In  [44]  a  total 
variational  type  optimization  problem  is  proposed  in  which  the  resulting  Euler-Lagrange 
equations  are  nonlinear  geometric  heat  equations  which  preserve  edges  much  better. 

The  optical  flow  constraint  above  is  of  course  very  strong.  Motivated  by  Moser  [51], 
we  have  proposed  a  modification  of  this  that  also  could  be  placed  in  a  variational  setting. 
Namely,  the  Moser  construction  described  in  [51]  allows  one  to  do  the  following:  Given 
a  family  of  nowhere-zero  2-forms  rt,  we  have  an  explicit  method  to  determine  a  family  of 
diffeomorphisms  cj)t  such  that 

$rt  =  r0. 

Differentiating  <f>lrt  =  r0  with  respect  to  t  yields 

Q 

Qjjt  +  (Vrt,  ut)  +  rtdiv(ut)  =  0.  (37) 

This  is  very  similar  in  form  to  the  standard  optical  flow  constraint  with  the  divergence 
term  added.  In  our  work,  we  interpret  image  intensity  as  a  type  of  form,  and  apply  the 
Moser  analysis  under  a  much  less  restrictive  assumption  than  the  standard  optical  flow 
constraint  given  in  equation  (36). 

6  Image  Interpolation  and  Optical  Flow 

We  would  like  to  show  how  optimal  transport  ideas  can  be  used  to  interpolate  a  given 
scalar  field  (and  even  a  vector  or  tensor  field)  from  one  image  to  another,  as  well  as  in 
optical  flow.  For  simplicity,  we  just  consider  the  intensity  maps  in  our  discussion  in  this 
section.  We  have  already  alluded  to  using  area-preserving  maps  to  weaken  the  optical 
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flow  constraint.  In  the  present  approach,  we  use  the  optimal  transport  philosophy  more 
directly. 

The  idea  is  to  minimize  a  functional  of  the  following  form  over  area-preserving  maps: 


M  — 


o  u-Gf 


(38) 


Here  the  first  term  controls  the  “goodness  of  fit”  between  the  images  G  :  flo  — *■  R  and 
F  :  — >  R,  and  the  second  controls  the  “smoothness”  of  the  mapping  u.  Some  other 
smoothness  term  could  be  used,  of  course.  As  with  the  Monge-Kantorovich  flows,  we  can 
assume  u  =  uo  s-1  where  u  is  an  initial  mapping  from  Oo  to  fii,  and  s  is  a  family  of 
area-preserving  mappings  from  Qo  to  itself  parameterized  by  time.  Since  we  are  dealing 
with  images,  Q0  and  Qi  will  both  be  the  same  rectangular  domain.  For  this  application, 
one  can  use  the  identity  map  for  u,  so  that  u  =  s-1.  We  also  assume  that  s  at  time  zero  is 
the  identity  map. 

The  methodology  given  above  for  optimal  transport  can  be  used  in  an  analogous  manner 
to  derive  local  and  non-local  schemes  for  decreasing  the  functional  (38).  One  can  exploit 
this  type  of  approach  as  the  basis  of  scalar  field  mapping  and  optical  flow  algorithms. 
Finally,  other  “goodness  of  fit”  terms  may  be  included.  In  fact,  we  plan  on  using  mutual 
information  in  this  context  making  contact  with  the  work  in  [72]  for  registration. 
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tor,  and  Jincheol  Ha),  to  appear  in  IEEE  Trans,  on  Aerospace  and  Electronic  Sys¬ 
tems. 

18.  “Flattening  maps  for  the  visualization  of  multi-branched  vessels”  (with  S.  Haker  and 
L.  Zhu),  to  appear  in  IEEE  Trans.  Medical  Imaging. 
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19.  “Statistically  based  flow  for  image  segmentation”  (with  E.  Pichon  and  R.  Kikinis), 
Medical  Image  Analysis  8  (2004),  pp.  267-274. 

20.  “The  3D  structure  of  real  polymer  foams”  (with  M.  Montminy  and  C.  Macosko), 
Journal  of  Colloidal  and  Interfacial  Science  280  (2004),  pp.  202-211. 

21.  “Area-based  medial  axis  of  planar  curves,”  (with  S.  Betelu,  M.  Niethammer,  and  G. 
Sapiro),  Int.  Journal  Computer  Vision  GO  (2004),  203-224. 

22.  “Flux  driven  automatic  centerline  extraction”  (with  S.  Bouix  and  K.  Siddiqi),  to 
appear  in  Medical  Image  Analysis,  2004. 

23.  “Statistical  analysis  of  RNA  backbone”  (with  E.  Hershkowitz,  G.  Sapiro,  and  L. 
Williams),  submitted  for  publication  to  IEEE  Trans.  Computational  Biology  and 
Bioinformatics. 

24.  “Dynamic  geodesic  active  contours”  (with  S.  Angenent  and  M.  Niethammer),  sub¬ 
mitted  for  publication  to  IEEE  Trans.  Automatic  control. 

25.  “Despeckling  of  medical  ultrasound  images”  (with  O.  Michailovich),  submitted  for 
publication  to  IEEE  Trans.  Medical  Imaging. 

26.  “Roles  of  site-bound  magnesium  ions  in  folding  of  globular  RNAs”  (with  E  Tan- 
nenbaum,  E.  Herschkowitz,  S.  Howerton,  G.  Perng,  and  L.  Williams),  submitted  for 
publication  in  RNA. 

27.  “Dynamic  denoising  of  tracking  sequences  of  SAR  data”  (with  O.  Michailovich), 
submitted  for  publication  to  IEEE  Trans.  Image  Processing. 

28.  “Fast  approximation  of  smooth  Functions  from  samples  of  partial  derivatives  with 
applications  to  phase  unwrapping”  (with  O.  Michailovich),  submitted  for  publication 
to  SIAM  J.  Numerical  Analysis. 

29.  “A  method  for  prediction  and  estimation  of  large-amplitude  optical  flows  -  an  ex¬ 
tended  Kalman  filtering  approach”  (with  O.  Michailovich),  submitted  for  publication 
to  IEEE  Trans.  Image  Processing. 

30.  “On  the  evolution  of  closed  curves  by  means  of  vector  distance  functions”  (with 
M.  Niethammer  and  P.  Vela),  submitted  for  publication  to  Int.  Journal  Computer 
Vision. 

31.  “Detecting  simple  points  in  higher  dimensions”  (with  M.  Niethammer,  W.  Kalies, 
and  K.  Mischaikow),  submitted  for  publication  in  Pattern  Recognition. 

32.  “Development  and  testing  of  highly  autonomous  unmanned  aerial  vehicles”  (with  J. 
Ha,  E.  Johnson,  and  A.  Proctor),  submitted  for  publication  in  AIAA  JACIC. 
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33.  “Automatic  recognition  of  electronic  components”  (with  G.  Pryor  and  A.  Goldstein), 
Circuits  Assembly  15  (2004),  pp.  36-40. 

34.  “Advanced  nonlinear  registration  algorithms  for  image  fusion”  (with  S.  Warfield  et 
at),  in  Brain  Mapping:  The  Methods,  Second  Edition  edited  by  Arthur  Toga  and 
John  Mazziotta,  Academic  Press,  pages  661-690,  2003. 

35.  “Mean  curvature  flows,  edge  detection,  and  medical  image  segmentation”  (with  S. 
Angenent,  S.  Haker,  and  A.  Yezzi),  in  Computational  Methods  in  Biophysics,  Bio¬ 
materials,  Biotechnology  and  Medical  Systems  edited  by  C.  Leondes,  pages  253-269, 
Kluwer,  2003. 

■  ' 

36.  “New  approach  to  Monge-Kantorovich  with  applications  to  computer  vision  and  im¬ 
age  processing”  (with  S.  Haker),  IMA  Series  on  Applied  Mathematics,  volume  133, 
Springer- Verlag,  New  York,  2003. 

37.  “Maximal  entropy  reconstruction  of  back  projection  images”  (with  T.  Georgiou  and 
P.  Olver),  IMA  Series  on  Applied  Mathematics,  volume  133,  Springer- Verlag,  New 
York  2003. 

38.  “Segmentation  of  diffusion  tensor  images”  (with  Eric  Pichon  and  Guillermo  Sapiro), 
in  Directions  in  Mathematical  Systems  Theory  and  Optimization  edited  by  Anders 
Rantzer  and  Chris  Byrnes,  pages  240-249,  Springer,  New  York,  2002. 

39.  “Optimal  image  interpolation  and  optical  flow”  (with  S.  Haker),  in  Multidisciplinary 
Research  in  Control,  Lecture  Notes  in  Control  and  Inform.  Sci  289,  pages  133-143, 
2002. 

40.  “Stochastic  crystalline  flows”  (with  G.  Ben-Arous  and  Ofer  Zeitouni),  in  Mathemat¬ 
ical  Systems  Theory  in  Biology,  Communications,  Computation,  and  Finance  edited 
by  J.  Rosenthal  and  D.  Gilliam,  IMA  Volumes  in  Mathematics  and  Its  Applications, 
volume  134,  pages  41-63,  Springer,  New  York,  2003. 

41.  “On  a  stochastic  model  of  geometric  snakes”  (with  D.  Nain,  G.  Unal,  A.  Yezzi,  and 
O  Zeitouni),  to  appear  in  Mathematical  Methods  in  Computer  Vision:  A  Handbook, 
edited  by  O.  Faugeras  and  N.  Paragios,  Springer- Verlag,  2005. 

42.  “Robust  control  and  tracking”,  Proceedings  of  IEEE  CDC’00. 

43.  “Affine  invariant  symmetry  sets”  (with  S.  Betelu  and  G.  Sapiro),  Proceedings  of 
ECCV’00 ,  Dublin,  Ireland,  June  2000. 

44.  “Nondistorting  maps  for  virtual  colonoscopy”  (with  S.  Angenent,  S.  Haker,  and  R. 
Kikinis),  Proceedings  of  SPIE,  San  Diego,  February  2000. 

45.  “New  approach  for  visualization  of  3D  colon  imagery”  (with  S.  Angenent,  S.  Haker, 
and  R.  Kikinis),  MICCAI’00,  October  2000. 
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46.  “Length-based  attacks  for  certain  group  based  encryption  rewriting  systems”  (with  J. 
Hughes),  SECI’02  Conference,  Securite  de  la  Communication  sur  Internet,  Septem¬ 
ber  2002. 

47.  “New  algorithms  for  3D  analysis  of  open-celled  foams,”  (with  M.  Montminy  and  C. 
Macosko),  Proceedings  of  FOAM  2000,  New  Jersey. 

48.  “High  resolution  sensing  and  anisotropic  segmentation  for  SAR  imagery”  (with  T. 
Georgiou),  Proceedings  of  IEEE  CDC’OO. 

49.  “Affine  invariant  erosion  for  3D  shapes”  (with  S.  Betelu  and  G.  Sapiro),  ICCV’01, 

2001. 

50.  “Missile  tracking  using  knowledge-based  adaptive  thresholding”  (with  S.  Haker,  G. 
Sapiro,  and  D.  Washburn),  ICIP’01,  2001. 

51.  “Cubical  homology  and  the  topological  classification  of  2D  and  3D  imagery”  (with 
M.  Allili  and  K.  Mischaikow),  ICIP’01,  2001. 

52.  “Optimal  transport  and  image  warping”  (with  S.  Haker),  IEEE  Conference  on  Vari¬ 
ational  and  Level  Set  Methods  in  Computer  Vision,  Vancouver,  20001. 

53.  “Mass-preserving  mappings  and  surface  registration”  (with  S.  Haker  and  R.  Kikinis), 
MICCAI’01,  October  2001. 

54.  “Minimal  transport  for  nonlinear  control”  (with  S.  Haker),  CDC’01,  December  2001. 

55.  “L1  based  optical  flow  for  cardiac  wall  motion  tracking”  (with  A.  Kumar,  S.  Haker, 
A.  Stillman,  C.  Curry,  D.  Giddens,  and  A.  Yezzi),  Proceedings  of  SPIE,  San  Diego, 
February  2001. 

56.  “Visual  tracking  and  object  recognition”  (with  A.  Yezzi  and  A.  Goldstein),  Proceed¬ 
ings  of  NICOLS’01,  St.  Petersburg,  Russia,  July,  2001. 

57.  “Conformal  flattening  maps  for  the  visualization  of  vessels”  (with  S.  Haker  and  L. 
Zhu),  Proceedings  of  SPIE,  San  Diego,  2002. 

58.  “Cubical  topological  analysis  of  blood  vessels”  (with  M.  Niethammer  and  A.  Stein), 
Proceedings  of  ICIP,  2002. 

59.  “Angle-reserving  mappings  and  multiply  branched  vessels”  (with  L.  Zhu  and  S. 
Haker),  Proceedings  of  ICIP,  2002. 

60.  “4D  active  surfaces  for  MR  cardiac  analysis”  (with  A.  Yezzi),  Proceedings  of  MIC- 
CAI’02. 

61.  “Flux  driven  fly-throughs,”  CVPR,  2003. 
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62.  “A  Stokes  flow  based  boundary  integral  formulation  for  measuring  crosssections  of 

two-dimensional  tubular  structures,”  (with  M.  Niethammer  and  E.  Pichonb  ICIP 
2003.  ' 

63.  “Algorithms  for  stochastic  approximations  of  curvature  flows”  (with  G.  Ben-Arous, 
N.  Shimkin,  G.  Unal,  and  O.  Zeitouni),  ICIP ,  2003. 

64.  “Dynamic  level  sets  for  visual  tracking”  (with  M.  Niethammer),  IEEE  Conference 
on  Decision  and  Control,  2003. 

65.  Statistically  based  surface  evolution  method  for  medical  image  segmentation:  pre¬ 
sentation  and  validation”  (with  E.  Pichon  and  R.  Kikinis),  MICCAI,  2003.  (Based 
student  presentation  award.) 

66.  “Active  contours  and  optical  flow  for  automatic  tracking  of  flying  vehicles”  (with  J. 
Ha,  C.  Alvino,  G.  Pryor,  E.  Johnson),  American  Control  Conference,  2004. 

67.  ’’Image  interpolation  based  on  optimal  mass  preserving  maps”  (with  L.  Zhu),  Pro- 
ceeedings  of  ISBI,  2004. 

68.  “Flux  driven  fly-throughs”  (with  S.  Bouix  and  K.  Siddiqi)  Proceedings  of  CVPR 
2003. 

69.  “Dynamic  geodesic  snakes”  (with  M.  Niethammer),  in  Proceedings  of  CVPR,  2004. 

70.  “Knowledge-based  3D  segmentation  and  reconstruction  of  left  coronary  arteries  using 
CT  images”  (with  D.  Giddens  and  Y.  Yang),  EMBSOf,  2004. 

71.  “Flow  patterns  and  wall  shear  stress  distributions  at  atherosclerotic-prone  sites  in 
a  human  left  coronary  artery-an  exploration  using  combined  methods  of  CT  and 
computational,  fluid  dynamics  (with  S.  Jin,  Y.  Yang,  J.  Oshinski,  A.  Tannenbaum 
J.  Gruden,  and  D.  Giddens),  EMBS04 ,  2004. 

72.  “Image  morphing  based  on  mutual  information  and  optimal  mass  transport”  (with 
L.  Zhu),  to  appear  in  Proceedings  of  ICIP,  2004. 

73.  “Automatic  tracking  of  flying  vehicles  using  geodesic  snakes  and  Kalman  filtering” 
(with  A.  Betser  and  P.  Vela),  IEEE  CDC,  2004. 

74.  “Flying  in  formation  using  a  pursuit  guidance  algorithm”  (with  A.  Betser,  G.  Pryor, 
and  P.  Vela),  submitted  for  publication  to  American  Control  Conference,  2005. 

75.  “Tracking  moving  and  deforming  shapes  using  a  particle  filter”  (with  Y.  Rathi,  N. 
Vaswani,  A.  Yezzi),  submitted  for  publication  to  CVPR,  2005. 

76.  “Affine  surface  evolution  for  3D  segmentation”  (with  Y.  Rathi,  P.  Olver,  G.  Sapiro), 
submitted  for  publication  to  ICIP,  2005. 
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77.  “Pattern  detection  and  image  segmentation  with  anisotropic  conformal  factors”  (with 
E.  Pichon),  submitted  for  publication  to  ICIP ,  2005. 


Books  Written  Under  AFOSR  Support 

1.  Invariance  and  Systems  Theory:  Algebraic  and  Geometric  Aspects,  Lecture  Notes  in 
Mathematics  845,  Springer- Verlag,  New  York,  1981. 

2.  Feedback  Control  Theory  (with  John  Doyle  and  Bruce  Francis),  MacMillan  Company, 
New  York,  1991.  (This  book  has  been  translated  into  Chinese  and  Japanese.) 

3.  Robust  Control  of  Distributed  Parameter  Systems  (with  Ciprian  Foias  and  Hitay 
Ozbay),  Lecture  Notes  in  Control  and  Information  Sciences  209,  Springer-Verlag, 
New  York,  1995. 

4.  Feedback  Control,  Uncertainty,  and  Complexity,  edited  by  Bruce  Francis  and  Allen 

Tannenbaum,  Lecture  Notes  in  Control  and  Information  Sciences  202,  Springer- 
Verlag,  New  York,  1995.  1 

5.  Deformation  Theory,  lectures  by  Michael  Artin,  notes  by  C.  S.  Seshadri  and  A. 
Tannenbaum,  Tata  Institute  Lecture  Notes,  Bombay,  India,  1976. 

6.  Curvature  Flows,  Visual  Tracking,  and  Computational  Vision ,  to  be  published  bv 
SIAM. 

7.  Mathematical  Methods  in  Computer  Vision,  edited  by  Peter  Olver  and  Allen  Tan¬ 
nenbaum,  IMA  Series  on  Applied  Mathematics,  volume  133,  Springer-Verlag,  2004. 


Patents  Based  on  AFOSR  Projects 

1.  “Conformal  Geometry  and  Texture  Mappings,”  (co-inventors  Sigurd  Angenent,  Steven 
Haker,  Allen  Tannenbaum,  and  Ron  Kikinis),  U.S.  Patent  Number  6,697,538,  issued 
February  24,  2004. 

2.  “4D  Kappa5  Gaussian  Noise  Reduction,”  (co-inventors  Harvey  Cline  and  Allen  Tan¬ 
nenbaum),  U.S.  Patent  Number  6,204,853  Bl,  issued  March  20,  2001. 

3.  “Curvature  Based  System  for  the  Segmentation  and  Analysis  of  Cardiac  Magnetic 
Resonance  Imagery,”  (co-inventors  Allen  Tannenbaum  and  Anthony  Yezzi),  U.S. 
Patent  Number  6,535,623  Bl,  issued  March  18,  2003. 
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1.  Allen  Tannenbaum  (faculty) 

2.  Marc  Niethammer  (Ph.D.  granted  December  2004) 

3.  Eric  Pichon  (Ph.D.  student;  expected  date  of  graduation:  September  2005) 

4.  Andrew  Stein  (M.S.  granted  December  2002) 
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